% Simulates data from the menu cost model with:
%
% 1. Constant aggregate demand
% 2. I.i.d. normal growth of the price level (non-zero mean)
% 3. AR(1)-normal idiosyncratic productivity shocks with homoskedastic errors
%
% Inputs:
%   NSim: Number of time periods to be simulated
%   burnin: Number of time periods at the start of the simulation that are
%      dropped
%   a: Grid for a (idiosyncratic tech shock)
%   rp: Grid for the real price
%   rho: Autoregressive coefficient in process for a_t
%   sigma_eps: Standard deviation of shocks to a_t process
%   mu: Average monthly inflation rate
%   sigma_eta: Standard deviation of shocks to the price level
%   F: The policy function
%   tfunc: If tfunc == 1 then the function prints out the time it takes to run
%   yesprint: If yesprint == 1 then the function prints 'progress reports'
%
% Outputs:
%   aSim: Simulated series for a (idiosyncratic tech shock)
%   InflSim: Simulated series for monthly inflation
%   rpSim: Simulated series for the real price
%   npSim: Simulated series for the nominal price
%   ncostSim: Simulated series for log(P_t)-log(A_it)
%   InflSim_sum: Simulated series for the nominal price level
%
% Notice: If this function starts printing out 'x' or 'y' it means that
% the grid for rp does not cover a large enough area.
%
% Jon Steinsson and Emi Nakamura, May 2006
%**********************************************************

function [aSim,InflSim,rpSim,npSim,ncostSim,InflSim_sum] ...
    = DataSim_hom_infl(NSim,burnin,a,rp,rho,sigma_eps,mu,sigma_eta,F,tfunc,yesprint)

time1 = cputime;

amax = a(end,1);
amin = a(1,1);
agridsize = size(a,1);

rpmax = rp(end,1);
rpmin = rp(1,1);
rpgridsize = size(rp,1);

if yesprint == 1
    fprintf('\n')
    fprintf('0 ')
end

% Simulate a time series for a
%******************************

aSim = zeros(NSim+burnin+1,1);
a_inc = (amax-amin)/(agridsize-1);
eps_a = sigma_eps*randn(NSim+burnin,1);
aSim(1,1) = a((agridsize/2 - mod(agridsize/2,1) + 1),1);
for i = 1:(NSim+burnin)
    delta_a = (rho-1)*aSim(i,1) + eps_a(i,1);
    delta_a_mod = mod(delta_a,a_inc);
    delta_a_mod_big = a_inc/2 < delta_a_mod;
    delta_a_mod_small = a_inc/2 >= delta_a_mod;
    delta_a_grid = delta_a - delta_a_mod_small*delta_a_mod ...
        + delta_a_mod_big*(a_inc-delta_a_mod);
    aSim(i+1,1) = aSim(i,1) + delta_a_grid;
    aSim(i+1,1) = min(amax,aSim(i+1,1));
    aSim(i+1,1) = max(amin,aSim(i+1,1));
end
aSim_num = int32((aSim - amin*ones(NSim+burnin+1,1))/a_inc+1);

if yesprint == 1
    fprintf('1 ')
 end
    
% Simulate a time series for Infl
%******************************

Infl_inc = (rpmax-rpmin)/(rpgridsize-1);
nu = sigma_eta*randn(NSim+burnin,1);
nu = [0; nu];
InflSim = mu + nu;
InflSim_mod = mod(InflSim,Infl_inc);
InflSim_mod_big = Infl_inc/2 < InflSim_mod;
InflSim_mod_small = Infl_inc/2 >= InflSim_mod;
InflSim = InflSim - InflSim_mod_small.*InflSim_mod ...
    + InflSim_mod_big.*(Infl_inc-InflSim_mod);
InflSim_num = int32(InflSim/Infl_inc);

if yesprint == 1
  fprintf('2 ')
end

% Simulate a time series for rp
%*****************************

rpSim = zeros(NSim+burnin+1,1);
rpSim_num = zeros(NSim+burnin+1,1);
rp_inc = (rpmax-rpmin)/(rpgridsize-1);
rpSim(1,1) = rp((rpgridsize/2 - mod(rpgridsize/2,1) + 1),1);
rpSim_num(1,1) = (rpSim(1,1) - rpmin)/rp_inc+1;
for i = 1:(NSim+burnin)
    % Adjust for inflation
    pre_adj_rp_num = rpSim_num(i,1) - double(InflSim_num(i+1,1));
    if pre_adj_rp_num > rpgridsize
        cutoff = mod(pre_adj_rp_num,rpgridsize);
        pre_adj_rp_num = rpgridsize;
        InflSim(i+1,1) = InflSim(i+1,1) - cutoff*Infl_inc;
        InflSim_num(i+1,1) = InflSim_num(i+1,1) - cutoff;
        fprintf('x ')
    end
    if pre_adj_rp_num < 1
        addon = 1-pre_adj_rp_num;
        pre_adj_rp_num = 1;
        InflSim(i+1,1) = InflSim(i+1,1) + addon*Infl_inc;
        InflSim_num(i+1,1) = InflSim_num(i+1,1) + addon;
        fprintf('y ')
    end
    
    % Calculate new price
    rpSim(i+1,1) = F(int32(pre_adj_rp_num),int32(aSim_num(i+1,1)));
    rpSim_num(i+1,1) = int32((rpSim(i+1,1) - rpmin)/rp_inc+1);
end

if yesprint == 1
  fprintf('3 ')
end

% Create Nominal Price Data
%**************************

npSim = zeros(NSim+burnin+1,1);
ncostSim = zeros(NSim+burnin+1,1);
InflSim_sum = zeros(NSim+burnin+1,1);
npSim(1,1) = rp((rpgridsize/2 - mod(rpgridsize/2,1) + 1),1);
ncostSim(1,1) = -a((agridsize/2 - mod(agridsize/2,1) + 1),1);
for i = 1:(NSim+burnin)
    InflSim_sum(i+1,1) = InflSim_sum(i,1) + InflSim(i+1,1);
end
npSim(2:end,1) = rpSim(2:end,1) + InflSim_sum(2:end,1);
ncostSim(2:end,1) = -aSim(2:end,1) + InflSim_sum(2:end,1);

if yesprint == 1
  fprintf('4 ')
  fprintf('\n\n')
end

if tfunc == 1
    fprintf('Time of DataSim: %8.3f\n\n',cputime - time1)
end

aSim = aSim((burnin+2):end,1);
InflSim = InflSim((burnin+2):end,1);
rpSim = rpSim((burnin+2):end,1);
npSim = npSim((burnin+2):end,1);
ncostSim = ncostSim((burnin+2):end,1);
InflSim_sum = InflSim_sum((burnin+2):end,1);
